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Unwanted motion of the probe with respect to the sample is a ubiquitous problem in 
scanning probe microscopy, causing both linear and nonlinear artifacts in experimental im¬ 
ages. We have designed a procedure to correct these artifacts by using orthogonal scan pairs 
to align each measurement line-by-line along the slow scan direction. We demonstrate the 
accuracy of our algorithm on both synthetic and experimental data and provide an imple¬ 
mentation of our method. 


Introduction 

Scanning probe microscopy (SPM) is a very pow¬ 
erful experimental tool due in part to the small 
probe size and corresponding high spatial reso¬ 
lution of contemporary physical and focused ra¬ 
diation probes. SPM experiments, such as scan¬ 
ning transmission electron microscopy (STEM), 
requires that the probe be moved across the sam¬ 
ple surface in order to form an image. This pro¬ 
cedure can introduce artifacts in the measure¬ 
ment due to the time delay between measure¬ 
ments and the accumulation of error in the probe 
position, for example from drift of the sample [1- 
8] . Examples of these artifacts include linear dis¬ 
tortions such as shear, expansion or contraction 
applied to the whole image, random “jitter” of 
each scanline’s origin position with respect to the 
intended position on the sample, jump disconti¬ 
nuities due to large sample jumps, and others. 
Virtually all scanning probe experiments contain 
image distortions, and these distortions are often 
large compared to atomic-scale features. 

Previous studies by other groups have attempted 
to measure and correct SPM distortions. Berkels 
et al. recorded and aligned series of STEM im¬ 
age exposures to improve peak precision [3] , and 
have also devised a method for nonlinear reg¬ 
istration of images series [5]. Jones and Nel- 
list corrected both linear and some nonlinear 
distortions in single STEM images by assum¬ 
ing prior knowledge of atomic features [4]. Sang 


and LeBeau have developed the “REV-STEM” 
method, where linear drift coefficients are mea¬ 
sured and corrected by the a series of STEM im¬ 
ages recorded in different scan directions [6]. In 
STEM experiments, many samples will be modi¬ 
fied or damaged by excessive electron dose [9, 10] 
and therefore using as few measurements as pos¬ 
sible to correct drift distortions is desirable. To 
that end, we have created a general case algo¬ 
rithm which requires no a priori assumptions 
about image features and requires only two im¬ 
ages as an input, minimizing acquisition time 
and electron dose needed for correction of drift 
artifacts. 

In this manuscript, we first show how nonlinear 
drift can affect scanning probe images. We then 
develop an algorithm to correct all linear and 
nonlinear drift distortions in scanning probe im¬ 
ages by correcting the scanline origin positions 
from two or more SPM images. We generate cor¬ 
rected images by using kernel density estimation 
to resample the images, and a develop a Fourier 
weighting scheme to further reduce error. We 
test this algorithm on synthetic data, and both 
simple and complex experimental datasets. We 
evaluate the algorithm by measuring deviations 
of atomic sites in corrected images from the best- 
fit lattice positions, and by measuring complex 
lattice strain fields before and after correction. 
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Theory 

The directionality of scanning probe 
microscopy 

In a standard scanning probe experiment, data is 
recorded pixel-by-pixel sequentially. The result¬ 
ing array consists of a “fast” direction, consisting 
of the primary travel direction of the probe, and 
the “slow” direction, along which the probe is 
moved after each line is completed. Because ad¬ 
jacent pixels in the fast scan direction are mea¬ 
sured much closer together in time, the informa¬ 
tion transfer along this direction is more reliable 
than the slow direction. After each line is com¬ 
pleted, the probe must be quickly moved over 
a long distance, and repositioned directly below 
the preceding line, along the slow direction. In¬ 
formation transfer along this direction is less re¬ 
liable, because it is recorded much more slowly, 
and depends on the ability of the probe to be po¬ 
sitioned with perfect accuracy after a change in 
probe speed and direction. If the sample or the 
electronic set-points of the microscope drift dur¬ 
ing the experiment, this motion will affect the 
slow direction much more strongly than the fast 
direction. 

In this paper, in order to develop an algorithm to 
correct scan distortions, we make the following 
assumptions about the scanning probe experi¬ 
ments: 


1. The scan directions are accurately known 
in the microscope frame of reference. 

2. The probe translation steps along the fast 
scan directions have negligible error. 

3. The sample is unchanged during or be¬ 
tween different scans. 


Given these assumptions, we will show that or¬ 
thogonal scans can be used to correct both linear 
and nonlinear scan distortions. 


The effects of linear and nonlinear drift 
distortion on scanning probe microscopy 

Fig. lA shows a synthetic image consisting of 
two different hexagonal lattices with an epitax¬ 
ial relationship. One lattice represents a sim¬ 
ple matrix, and the other forms a moire pattern 
with the matrix, inside a simulated precipitate, 
forming an image with both low and high spa¬ 
tial frequencies. In an ideal scanning probe ex¬ 
periment, the sample remains perfectly station¬ 
ary during collection of the data, and the probe 
is perfectly positioned with respect to the sam¬ 
ple. An ideal resampled scanning probe image is 
shown in Fig. IB, with a horizontal fast scan di¬ 
rection. This image is a perfect representation of 
the input data, since the sampling rate is below 
the Nyquist limit. 

Scan distortions along the slow scan direction 
can be equivalently represented by either mov¬ 
ing the sample with respect to an ideal scan, 
or by moving the scanline origins with respect 
to the sample. Moving the sample continuously 
during a scanning probe experiment produces a 
linear drift artifact, like that shown in Fig. IC. 
This drift artifact can be described by a linear 
affine transformation of the underlying data [11], 
consisting of an expansion in one direction and a 
compression in another. Drift of this type is rou¬ 
tinely corrected in scanning probe experiments 
[6, 12, 13]. 

Examples of more complex drift artifacts are 
shown in Fig. ID. The artifacts could be caused 
by sudden jumps in the sample position, jumps 
in the probe positioning system, charging of the 
sample over time, or changes in the sample-probe 
environment. In this example, the drift direc¬ 
tion changes several times during data collection. 
Additionally, a jump discontinuity is present 
roughly one third from the top of the image. The 
simulated SPM image shows large distortions, 
and the inset Fourier transform shows streaks 
along the slow scan direction. The streaks are 
caused by the various locally-distorted regions 
of the sample being misaligned with respect to 
each other. 
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FIG. 1. (A) Synthetic dataset with examples of 

scanning probe microscopy images recorded with and 
without error in the probe positions, depicted by the 
colored outline where the scanline origins are repre¬ 
sented by a thicker line. The cases considered are 
(B) perfect sampling, (C) linear drift distortion, (D) 
nonlinear drift distortions and a jump discontinuity, 
and (E) random noise added to the scanline origins. 
The square root of the Fourier transform amplitude 
is inset into each image. 


Between scanlines, the probe must be reposi¬ 
tioned directly adjacent to the origin of the pre¬ 
vious scanline with high precision. If there is 
some error (perhaps due to electronic noise or 
hysteresis), it could manifest as a random devi¬ 


ation between the ideal and the actual scanline 
origins. An example of this is shown in Fig. IE, 
which is very similar to flagging often observed 
in scanning probe microscopies. This artifact 
also produces strong streaks along the slow scan 
direction in the Fourier transform. The combi¬ 
nation of the above artifacts can represent many 
nonlinear drift distortions observed in experi¬ 
ments. In the next sections, we will show that 
these distortions can be corrected with high ac¬ 
curacy. 


An algorithm to correct nonlinear distortions 
using orthogonal scan directions 

As in the previous section, we will represent all 
image distortions by translation of the scanline 
origins. In order to reverse these distortions, 
we therefore need to estimate all scanline origins 
from the distorted data. From a single scan, this 
is impossible without making some assumptions 
about the underlying data (such as assuming a 
periodic lattice). However, we can make use of 
the accurate information transfer along the fast 
scan direction, which we assume to be essentially 
error-free. By recording multiple images along 
orthogonal directions, we can use the fast scan 
direction of one image to calibrate the slow scan 
direction of another. This procedure is outlined 
in Fig. 2. 

This algorithm requires a minimum of two scan¬ 
ning probe images, ideally recorded with fast 
scan directions 90 degrees apart. We begin by 
initializing a reconstruction space that is large 
enough to contain all transformed and translated 
images. Samples with a large amount of drift 
will require a larger padding region. The images 
are first rotated to the specified scan directions, 
and roughly aligned using phase correlation. Af¬ 
ter applying a low pass filter, a reference point 
is selected where the mean absolute difference 
between all images reaches a minimum. 

The next step is a rough initial alignment of all 
images. Each image is aligned relative to a strip 
of user-defined width taken from whichever other 
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■ Import all image datasets, scan directions 0, 
and step sizes along the scan lines. 
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■ Generate padded image space, and initialize 
all image scanline origins and directions. 

■ Use phase correlation or cross correlation 
for initial alignment guess. 

■ Determine initial reference location as the 
point with the smallest absolute difference. 


For each image, select a reference strip from 
other image at most orthogonal direction. 

Use linear fit of fitted origins to calculate 
initial step vectors to adjacent origins. 

Along successive lines, update origin position 
by testing it on adjacent pixels, using total 
absolute difference between scanline 
measurement and interpolation of reference. 

Loop through all scan lines in all images, 
a user-specified number of times. 


Use phase correlation or cross correlation 
to check global alignment between images. 


For each image, calculate reference image 
using bilinear interpolated images calculated 
from current origins for all other images. 

Along successive lines, update origin position 
by testing steps in cardinal directions, using 
total absolute difference between scanline 
measurement and interpolation of reference. 



■ If scanline origin does not move, 
halve the search step size. 


yyes 


Loop through all scan lines in all images, 
and then either perform global 
alignment, or terminate reconstruction. 


Generate corrected images for each scan 
direction using kernel density estimation on 
an upsampled image space. 

Combine images in Fourier space using 
weighting cos2(0), and masking in real space. 


FIG. 2. Schematic of the nonlinear correction algo¬ 
rithm. 


image was recorded at a direction closest to 90 
degrees away, starting from the reference point 
and moving outwards. Each scanline is scored 
by taking the mean absolute difference between 
the measured intensity and a resampled estimate 
of the reference strip. Next the scanline origin is 
moved to the four adjacent pixels and the associ¬ 
ated score is computed. Whichever scanline ori¬ 
gin produces the lowest score is selected, and the 
origin position is updated. This update is stabi¬ 
lized by under-relaxation, i.e. rather than mov¬ 
ing one pixel in a given direction, the scanline 
will be moved in that direction by some reduced 
factor. As the alignment proceeds outwards from 
the reference point, new scanline origins are gen¬ 


erated by linearly fitting the step size between all 
fitted origins, and adding this step to adjacent 
origins sequentially. The purpose of this rough 
alignment is to robustly estimate the global lin¬ 
ear drift between different images and provide 
a good starting guess for the next step. These 
iterations are repeated several times and can be 
checked visually by the user at this stage. 

After the initial alignment is complete, we it¬ 
eratively perform the main alignment of all 
datasets. First, a synthetic image is generated 
for each recorded dataset using kernel density es¬ 
timation (KDE), which is described in the Ap¬ 
pendix. After a synthetic image has been gen¬ 
erated from all datasets, the primary refinement 
proceeds. Each dataset is aligned to the average 
of all other datasets (at least one, but could be 
any number), one scanline at a time. This refine¬ 
ment is similar to that described above, where 
each scanline origin is scored using mean abso¬ 
lute difference and compared to an origin shift in 
each of the 4 nearest-neighbor pixel directions, 
with the shift size typically initialized to 0.5 pix¬ 
els. However in this refinement step, if a given 
scanline origin does not move, its search step size 
is halved. After all scanlines have been refined 
once, the synthetic images are recalculated. At 
this stage, global alignment between all images 
can be checked via phase correlation. If an im¬ 
age shift of more than a few pixels is detected, 
the origin shift step size is reset to the initial 
value for that image. These steps are repeated 
until either the total origin shift falls below a 
convergence threshold (typically 1 pixel per im¬ 
age) , or a user-specified number of iterations are 
completed (~20 iterations required for a typical 
experimental STEM dataset). 

At this stage the user again reviews the trans¬ 
formed images, and if required performs more 
primary refinement iterations. Once the refine¬ 
ment is completed, the final step is to gener¬ 
ate an output image. This can be done simply 
by taking the mean of the transformed images. 
However, this method is not ideal; the nonlin¬ 
ear drift correction algorithm described above 
cannot fully reconstruct the original data since 
some information can be irrevocably lost along 
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the slow scan direction, due to noise or large 
nonlinear drift distortions. A better method is 
to weight the fast scan direction of each trans¬ 
formed image more strongly than its slow scan 
direction. This is accomplished by taking the 
Fourier transform of each image Jj, and multi¬ 
plying by the weighting function cos^ 9 where 
9 is the angle between each Fourier coordinate 
pixel and the scan direction. The N weighted 
images are then added together, and divided by 
the sum of all weighting functions. The final im¬ 
age /output is given by the inverse Fourier trans¬ 
form of the weighted sum: 


^|/< 




r_' 


( 1 ) 


where ^{} represents the 2D Fourier transform. 
This weighting can be tightened further towards 
the fast scan direction if more scan angles are 
recorded between 0 and 90 degrees. This weight¬ 
ing can improve the corrected image quality by 
suppressing noise in the slow scan direction. 


Our method has several limitations: perform¬ 
ing drift correction line-by-line requires that the 
data contains enough variance along each line 
to be accurately aligned. If the data does not 
contain sufficient high spatial frequency infor¬ 
mation, additional constraints must be applied. 
One potential constraint would be to apply a 
moving average to the scanline origins, or to 
force the origins to follow a linear path. We 
therefore recommend images at atomic resolu¬ 
tion be recorded at a slight angle with respect 
to crystallographic planes. Another potential 
limitation occurs when applying our method to 
nearly-perfect periodic lattices, such as atomic- 
resolution images aligned to a low-index crys¬ 
tallographic zone axis. A global lattice mis¬ 
alignment between images can easily occur, and 
therefore we require enough low spatial fre¬ 
quency information (such as an edge, a precipi¬ 
tate or alignment markers) to detect and correct 
lattice misalignments. High-dose STEM experi¬ 
ments can also erode or damage the sample be¬ 
ing imaged [9, 10]. Therefore the electron dose 
should be sufficiently low for the sample to re¬ 
main unchanged between scans. Finally, we have 
assumed near-perfect information transfer along 


the fast scan direction of each scanline. Sys¬ 
tematic errors along this direction, for example 
electronic noise at 60 Hz in the STEM scanning 
coils, will cause misalignment of the relative im¬ 
ages. 


Implementation of our algorithm in 
MATLAB 


We have provided an implementation of our al¬ 
gorithm programmed in MATLAB in the sup¬ 
plementary materials of this manuscript. This 
implementation should be general enough to cor¬ 
rect any experimental datasets, and contains 
many user-selectable settings to aid in recon¬ 
structions. The code has been tested on datasets 
consisting of up to 500 images. It has also 
been tested for simultaneously recorded images, 
for example bright field and dark field pairs in 
STEM. A synthetic nonlinear drift example like 
that described below is also included. 


Results and Discussion 


Synthetic dataset example 


We first test the above algorithm using simu¬ 
lated image pairs. Simulated tests are important 
because unlike in a real experiment, we know 
all of the original data values and can therefore 
precisely assess any residual error after correc¬ 
tion. The synthetic image data is plotted in 
Eig. 3A, with the scanline origins for an im¬ 
age pair with a large amount of nonlinear drift 
overlaid. The progression of our drift correction 
algorithm applied to this dataset is shown in 
Eig. 3B. After the initial alignment step to es¬ 
timate image translation, several different moire 
lattice patterns are visible on the left image of 
Eig. 3B. The rough alignment step removes the 
long-range nonlinear distortions, producing the 
smooth outer lattice shown in the center of Eig. 
3B, though with a unit cell jump between the 
two images (visible in the precipitate in the cen- 
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FIG. 3. (A) Synthetic image with scanline origins containing nonlinear drift shown as a colored outline, for a 
horizontal scan direction (red) and a vertical scan (green). (B) The three stages of drift correction applied to 
the example in (A). Drift correction examples for the data in (A) are shown with (C) high S/N (D) low S/N, 
and (E) scans with random Gaussian shifts (2 pixel standard deviation) of each scanline origin (“jitter”) 
added to the scanline origins shown in (A). The first three rows of (G)-(E) show original and reconstructed 
images, bottom row shows residual noise after subtracting original signal with a range of ±25% normalized 
intensity. Windowed Fourier transforms inset into all images. 


ter). Finally, the main alignment step corrects 
the unit cell jump and removes all residual drift. 
This can be seen by comparing the scanline ori¬ 
gins in Fig. 3B to those of the original samples 
in Fig. 3A. 

Two outputs of our nonlinear drift correction al¬ 
gorithm are shown in Figs. 3C-D, corresponding 
to high signal-to-noise (S/N) images, and low 
S/N images respectively. In both cases, the re¬ 
constructed image is very accurate and contains 
no obvious nonlinear drift artifacts. The high 
S/N intensity residual is uniformly random and 
contains almost no residual signal of the crystal 
lattices. The low S/N reconstruction has a sim¬ 
ilarly uniformly random residual. However, the 
inset Fourier transform in Fig. 3D shows that 
the noise signal is biased along more horizon¬ 
tal and vertical directions. This is because the 


noise is more correlated along the image pair 
scan directions, meaning that any error in the 
corrected scanline position will couple along the 
entire scanline. Despite this issue, the recon¬ 
structed image is still a very good match to the 
original data, shown by the lack of residual crys¬ 
talline peaks in the residual Fourier transform. 

The drift correction example shown in Fig. 3E 
has infinite S/N, but contains random shifts of 
all scanline origins, with a magnitude approx¬ 
imately half the atomic spacing. The recon¬ 
structed image closely matches the original syn¬ 
thetic data. However, the residual image shows 
that the signal contrast (peak-to-valley ratio) is 
decreased, relative to the original data. This is 
caused by the interpolation step of the recon¬ 
struction algorithm applied to the highly non¬ 
linear sampling of the synthetic data. Interpo- 
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lation tends to reduce the deviation of an image 
from the local mean, leading to a residual signal 
in Fig. 3E approximately equal to the intensity 
mean minus the original data multiplied by a 
small negative value. Despite this, the recon¬ 
structed image is still a good match to the orig¬ 
inal data, with the peak locations unchanged. 


Simple experimental example 

We next test our algorithm using a simple ex¬ 
perimental dataset, where the sample consists of 
two perovskite layers, a Lao. 2 Bio. 8 Fe 03 film on a 
SrRuOa buffer electrode, imaged along the [001] 
zone axis. The detector used creates a high angle 
annular dark field (HAADF) image that is sensi¬ 
tive to atomic number. The Lao. 2 Bio. 8 Fe 03 layer 
appears brighter due to the high-Z elements La 
{Z = 57) and Bi (Z = 83) on the perovskite 
A-site. The SrRu 03 layer appears dimmer, as 
its brightest peak occurs for Ru (Z = 44) on 
the perovskite B-site. Two HAADF images were 
recorded with orthogonal scan directions, and 
the resulting corrected image is plotted in Fig. 
4A. The Fourier transform amplitudes for the 
original vertical scan image, the original hori¬ 
zontal scan image and the corrected image are 
plotted in Figs. 4B, E, and H respectively. As 
in the synthetic example of the previous section, 
the correction routine has removed the streaks 
along the slow scan directions. 

To quantify the improvement in data quality af¬ 
ter reconstruction, we have performed nonlinear 
least squares fits of 2D Gaussian peaks (5 param¬ 
eter fits: peak position, peak intensity, standard 
deviation and a constant background intensity) 
to all atomic columns in a 50 x 50 unit cell re¬ 
gion, shown in Fig. 4A. We compared all column 
position fits to a global best-fit linear lattice, and 
calculated the deviation ctrms from the ideal po¬ 
sitions. A 2D histogram of the La/Bi and Fe 
sites for the vertical scan image is plotted in Fig. 
4C and D, for the horizontal scan image in Fig. 
4F and G, and for the corrected image in Fig. 41 
and J respectively. 


The results show that for both original images, 
long-range nonlinear distortions have produced 
large asymmetric deviations from the ideal lat¬ 
tice. By comparison, the corrected image devia¬ 
tions form a tighter and more symmetric pattern 
of deviations from the ideal lattice. The root- 
mean-square of all deviations <trms for the La/Bi 
and Fe sites has been reduced by a factor of ap¬ 
proximately three and two respectively. More 
importantly, the asymmetric deviations caused 
by sample drift distortions have been almost 
completely removed. 


Complex experimental example 

In the previous section, the sample consisted of 
a nearly ideal lattice, which could have been ap¬ 
proximately corrected by assuming an ideal lat¬ 
tice. To demonstrate that our algorithm works 
with arbitrarily complex experimental data, we 
now examine an experimental dataset contain¬ 
ing a large amount of nonlinear distortions over 
many different length scales. This sample is 
composed of a superlattice of alternating per¬ 
ovskite layers, similar to those described by 
Schlom et al. [14], with a non-uniform lattice 
structure as shown in Fig. 5 and detailed in 
Ref. [15]. This sample was chosen because these 
HAADF images contain a large amount of non¬ 
linear drift distortion, and it is known to contain 
a complex strain field. Again, we have recorded 
orthogonal scan pairs and applied our algorithm 
to correct the drift distortions, using directional 
Fourier weighting to produce the final image. 
Strips from the same region of the original im¬ 
ages and the corrected image are plotted in Figs. 
5A, B and C respectively. A jump discontinuity 
is visible running across the left side of Fig. 5B, 
but it has been removed from Fig. 5G. Both of 
Figs. 5A and B contain small long-range lattice 
distortions that are not present in Fig. 5C. These 
distortions could easily be mistaken for real vari¬ 
ations in the sample without application of non¬ 
linear drift correction. 

We have performed a measurement of the lat¬ 
tice strains from both the A (brighter) and the 
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FIG. 4. (A) Drift-corrected HAADF image of a Lao, 2 Bio, 8 Fe 03 film from an orthogonal scan pair, imaged 
along the [001] direction. Fourier transform amplitudes of (B) vertical fast scan direction image, (E) hor¬ 
izontal fast scan direction image, and (H) drift-corrected image given in (A). Two-dimensional deviations 
from best-fit lattice for box shown in (A) for the La/Bi sites are plotted in (C), (F), and (I), and deviations 
of the Fe sites in (D), (G), and (J), in the same order as (B), (E) and (H). RMS 2D deviation ctrms shown 
for all sites/images. 


B (dimmer) perovskite sites, for both of the orig¬ 
inal images and the corrected image. The lattice 
strains are calculated using a real space fitting 
procedure. First, all site positions are measured 
using a 2D Gaussian fit. Then, a best-fit ref¬ 
erence lattice was calculated for all sites, and 
the deviations from the ideal lattice positions 
are computed. These deviations are transformed 
into a 2D image using KDE, as in Fig. 6, with 
a kernel standard deviation of 1 unit cell length. 
Finally, the 2D strain maps are calculated by nu¬ 
merical differentiation of the displacement maps. 
The site intensities and strain maps for the orig¬ 
inal images and the corrected image are shown 
in Figs. 5D-F. 

The strain maps in Figs. 5D and E show clear 
signs of drift distortion, manifesting as long 
streaks along the slow scan directions that span 
the field of view. These streaks are espe¬ 
cially pronounced in strain measurements, since 
numerical differentiation magnifies any noise 
present. The only strain maps that do not show 
large streak artifacts are those measured along 
a direction close to the fast scan direction, Cxx 


for Fig. 5D, and Cyy for Fig. 5E. By contrast, no 
artifacts are visible in the drift-corrected strain 
maps plotted in Fig. 5F. This example demon¬ 
strates that our algorithm can successfully cor¬ 
rect nonlinear drift without removing complex 
real signals in the underlying data. A detailed 
analysis of these signals is given in Ref. [15]. 


Conclusion 

We have described a procedure for correcting 
nonlinear drift distortion in scanning probe im¬ 
ages using orthogonal scan pairs. Our algorithm 
corrects drift distortion in each image line-by- 
line, by moving each scanline origin to minimize 
the absolute difference between the measured 
line intensity and the current iteration of the cor¬ 
rected image recorded along an orthogonal direc¬ 
tion. We have successfully applied the algorithm 
to synthetic data, and both simple and complex 
experimental datasets. Our results show that 
recording only two images is sufficient to remove 
most of the nonlinear drift error. Including ad- 
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-4% 0% 4% -4% 0% 4% 


FIG. 5. Example of nonlinear drift correction for a complex experimental dataset. Strip of a square HAADF 
image for (A) near-vertical scan direction, (B) near-horizontal scan direction and (C) drift-corrected image 
from (A) and (B). (D), (E), and (E) rows contain Eourier transform amplitudes, site intensities and strain 
field measurements for A and B sites, measured from (A), (B) and (C) respectively. 


ditional images in the reconstruction will further 
reduce the error, especially from additional scan 
directions. Our algorithm should prove useful 
for all scanning probe experiments where high 
precision is required. 


Acknowledgements 

We thank Wolfgang Theis, Peter Ercius, Mary 
Scott and Matt Bowers for helpful discussions. 
Work at the Molecular Foundry was supported 
by the Office of Science, Office of Basic Energy 
Sciences, of the U.S. Department of Energy un¬ 
der Contract No. DE-AC02-05CH11231. 


Appendix 

Kernel Density Estimation 

Kernel density estimation (KDE) is a simple and 
robust method to estimate the values of a con¬ 
tinuous field from non-uniformly spaced discrete 


measurements [16, 17]. The procedure we use is 
demonstrated in Fig. 6. All measurement pix¬ 
els from all scanlines are first added to the re¬ 
construction space using bilinear weighting of a 
2x2 pixel region. Then, KDE is used to inter¬ 
polate the value of all image pixels. The “band¬ 
width” of the KDE procedure is set by the stan¬ 
dard deviation cr of a 2D Gaussian kernel, and 
is typically set to a value of 0.5 pixels for this 
algorithm. If cr is much less than the distance 
between samples, KDE will asymptotically ap¬ 
proach a nearest-neighbor interpolation. As cr is 
increased, the interpolation becomes smoother, 
shown in Fig. 6. This method is very robust 
against large nonlinear jumps between adjacent 
scanlines, where no measurement estimate may 
be available. Equally important is the resulting 
image and its derivatives are free of discontinu¬ 
ities, which is important for peak fitting or strain 
measurements. For comparison, a linear interpo¬ 
lation calculated on a nonuniform triangular grid 
is plotted the upper right plot of Fig. 6. This in¬ 
terpolation produces a poor result with a slope 
discontinuity at each triangular boundary. 
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FIG. 6. Top row shows (from left to right) list of 
peaks with locations r, bilinear weighting of each 
sample, bilinear weighting times sample intensity Ip, 
and a linear interpolation of these points using trian¬ 
gular regions. All other rows show grid interpolation 
of discrete samples using kernel density estimation 
(KDE), with a different kernel size a for each row. 
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